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ABSTRACT 

In a companion paper ( Laibe & Price|20lTb i, we have presented an algorithm for simulating 



two-fluid gas and dust mixtures in Smoothed Particle Hydrodynamics (SPH). In this paper, 
we develop an implicit timestepping method that preserves the exact conservation of the both 
linear and angular momentum in the underlying SPH algorithm, but unlike previous schemes, 
allows the iterations to converge to arbitrary accuracy and is suited to the treatment of non- 
linear drag regimes. The algorithm presented in Paper I is also extended to deal with realistic 
astrophysical drag regimes, including both linear and non-linear Epstein and Stokes drag. The 
scheme is benchmarked against the test suite presented in Paper I, including i) the analytic 
solutions of the dustybox problem and ii) solutions of the dustywave, dustyshock, dustysedov 
and DUSTYDisc obtained with explicit timestepping. We find that the implicit method is 1-10 
times faster than the explicit temporal integration when the ratio r between the the timestep 
and the drag stopping time is 1 < r < 1000. 

Key words: hydrodynamics — methods: numerical — ISM: dust, extinction — protoplane- 
tary discs — planets and satelhtes: formation 



1 INTRODUCTION 

Dust in cold astrophysical systems spans a huge range of sizes 
from sub-micron sized grains in the interstellar medium to kilo- 
metre sized planetesimals involved in planet formation. Moreover, 
the ratio of dust to gas as well as the density and temperature of 
the gaseous environment in which dust is embedded can also vary 
strongly. Handling this full range of physical parameters presents 
a challenge to numerical schemes designed to simulate dusty gas 
in astrophysics. The main challenges are i) that at high drag (e.g. 
small grains), the small timestep required means that purely explicit 
timestepping methods become prohibitive and ii) that a wide range 
of physical drag prescriptions, including non-linear drag regimes, 
need to be handled by the code. 

Two main prescriptions for drag between gas and solid parti- 
cles are applicable to astrophysics: The Epstein regime — where 
the gas surrounding a grain can be treated as a dilute medium — 
and the Stokes regime — where the grains can be treated as solid 
bodies surrounded by a fluid — (see e.g. |Baines et al.|1965||Stepin^ 
|ski & Valageas|1996 t. The dependence of the drag term on local 
parameters of the gas (density, temperature) and the dust (typical 
grain size, mass) differ between the two regimes, in turn implying 
very different dynamics for the dust grains. For example, in a proto- 
planetary disc, both of these regimes may be applicable in different 
regions of the disc. 

In a companion paper ( [Laibe & Price|[2bl lb[ hereafter Pa- 
per I), we have developed a new algorithm for treating two-fluid 
gas-dust astrophysical mixtures in Smoothed Particle Hydrody- 



namics (SPH). Benchmarking of the method demonstrated that the 
algorithm gives accurate solutions on a range of test problems rele- 
vant to astrophysics and substantially improve previous algorithms 
(Monaghan & Kochar yan| 1995| l. However, in Paper I, we used only 
a simple explicit time stepping and considered only linear drag 
regimes with a constant drag coefficient. In this paper, we present 
an implicit timestepping method that can be applied to both lin- 
ear and non-linear drag regimes, which is both more accurate and 
more general than the scheme proposed by Monaghan] ^1997^ . We 
also discuss the SPH implementation of both the Epstein and the 
Stokes regimes in their full generality. 



The paper is organised as follows: The equations of motion 
and the characteristics of the different astrophysical drag regimes 
for gas and dust mixtures are given in Sec. [2] We summarise the 
SPH formalism used for integrating these equations (derived in de- 
tail in Paper I) in Sec.[3]and extend it to deal with the drag regimes 
encountered in astrophysics. Particular attention is paid in Sec. [4] 
to improving the Monaghan ( 19971 implicit timestepping scheme, 
including its generalisation for non-linear drag regimes. Finally, 
the algorithm for non-linear drag regimes is tested against the ana- 
lytic solutions of the dustybox problem, as well as the dustywave, 
DUSTYSEDOV, DUSTYSHOCK and DUSTYDISC tests, in Sec.|5] 
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2 GAS AND DUST EVOLUTION IN ASTROPHYSICAL 
SYSTEMS 

2.1 Evolution equations 

The equations describing the evolution of astrophysical gas and 
dust mixtures, where dust is treated as a pressureless, inviscid, con- 
tinuous fluid have been described in detail in Paper I. The equations 
in the continuum limit are given by: 
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where the subscripts g and d refer to the gas and dust, respectively 
such that Pg is the gas pressure, Vg and vj are the fluid velocities and 
u is the specific internal energy of gas. The volume densities of gas 
and dust (pg and pd, respectively) are related to the corresponding 
intrinsic densities (pg and pd, respectively) according to 



Pd = (1 - d)pA, 
Pg = fPg' 
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where 6 is the volume filling fraction of the dust. Finally, the drag 
force and heating terms are given by: 
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The drag coeflicient K has dimensions of mass per unit volume per 
unit time and is generally a function of the relative velocity between 
the two fluids Av = |Vg- Vdl, implying a non-linear drag regime with 
respect to the differential velocity between the gas and the dust. In 
most of the astrophysical systems, the dust is diluted enough into 
the gas so that the gas filling fraction is 1 to a very good level 
of approximation, such that the dust buoyancy term (1 - d)VP, is 
negligible. 



2.2 Astrophysical drag regimes 

Microscopic collisions of gas molecules on a single dust grain re- 
sult in a net exchange of momentum which is equivalent to a drag 
force Fdrag between the two phases. Two limiting cases occur when 
comparing the typical geometrical size s of a dust grain to the mean 
free path A, of the gas. 

When the typical grain size is negligible compared to the col- 
lisional mean free path of the gas particles (s <c A„), the grains 
are surrounded by a dilute gas phase and may be treated using 
the Epstein drag prescription. In this limit, the analytic expression 
of the resulting drag force has been derived, assuming spherical, 
compact grains with homogeneous composition, for both specular 
and diffuse reflections on the grain surface (see |Baines et al.|1965| 
for the complete derivation). These expressions have been widely 
used in astrophysical studies (see e.g. [Chiang & Youdin|2010| for 
references), sometimes incorrectly where the grains are known to 



be porous and have fractal structures l |Blum & Wurm|2008[ l (thus 
breaking the assumptions of the Epstein prescription). 

For grain sizes larger than the collisional mean free path 
(s » /ig), grains experience a local differential velocity with re- 
spect to a uniform viscous flow and should be treated using the 
Stokes drag prescription ( [Fan & Zhu|I998| . In this case, the mo- 
mentum is diffused by viscosity into the fluid, which implies that 
the drag expression strongly depends on the local Reynolds number 
defined according to 



2s Vd - Vg 



(10) 



where v is the kinematic viscosity of the gas. Analytic expressions 
for the drag force can be derived at small Reynolds numbers. At 
higher Reynolds number, the drag law is inferred from experiments. 
Rigorously, additional contributions to the drag should arise from 
the grain acceleration (carried mass and Basset contribution), grain 
rotation (Magnus) and in the presence of strong local shear, pres- 
sure and temperature gradients (e.g. [Fan "& Zhu 1998 1. These cor- 
rections are negligible in nearly all astrophysical contexts. 

No current analytic theory describes how both the gas and 
the dust fluid exchange momentum in the intermediate regime (i.e. 
.s - /ig). Generally, an asymptotic continuous interpolation between 
the two limiting Epstein and Stokes regimes is used. [Stepinski &[ 



Valageas 1 1996 1 suggest adopting s 



i/ig as a means of obtaining 
a smooth transition. It should be noted that although this approach 
is convenient, there is no clear measure of the physical accuracy of 
this assumption. 

2.2.1 Epstein regime for dilute media 

In a dilute medium (A, > 4s/9), grains are small enough not to 
disturb the Maxwellian distribution of the gas velocity. Assuming 
grains are spherical, that the mass of a gas molecule is negligible 
compared to the mass of a dust grain and that the reflection of gas 
particles from collisions with dust grains are specular, the expres- 
sion of the drag force on a single grain Fdrag (which differs from the 
volume force F^^^^ by a factor Pi/m^, see Paper I) for the Epstein 
regime is given by 
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where s corresponds to the grain radius and m, p„, T denote the 
mass of the gas molecules, the intrinsic gas density and the local 
temperature of the mixture (the gas and the dust are supposed to 
have the same temperature). The thermal sound speed of the gas 
is thus Cs = -sJyk^T/m and the mean thermal velocity of the gas, 
t's -\/8M7- The dimensionless quantity if/ is defined according to 



y Av 
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where Av = Vj- Vg = Av x is the differential velocity (x being a unit 
vector). However, depending on the characteristics of the problem 
(i.e. low or high Mach numbers, or both), simpler and computa- 
tionally less expensive approximations may be used. For i/r <c I, 
i.e. low Mach numbers, Eq.[TT|can be expanded to third order in if/, 
giving 
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which is usually simplified to its linear term, 
An 



c'sAv. 
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For 1^ » 1, i.e. high Mach numbers, the Taylor expansion in 
of Eq.|l 1 [gives 
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which is usually reduced to its quadratic term. 
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A convenient way to handle Epstein drag at both low and high 
Mach numbers is to use an interpolation between the two asymp- 
totic regimes given by Eqs. 14 and|16| as derived in |Kwok|jl975^ 
(cf . [Paardekooper & Mellema|2006^ , giving 
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The deviation of Eq.[T7]from the full expression (Eq.[TTJ is < 1% 
( |Kwok|[T975[ >. Thus, in general, we adopt Eq. [TT] for the Epstein 
regime. We compare the differences between the various Epstein 
expressions in Sec.|5] 

2.2.2 Stokes regime for dense media 

In a dense medium (/!„ > 4i/9), grains should be treated with the 
Stokes drag regime, for which the expression of the drag force Fj^g 
is: 
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where the coefficient Co is a piecewise function of the local 
Reynolds number: 



Cd 
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where is defined in Eq. [TO] Equation [T9] indicates that at small 
Reynolds numbers {R^ < 1), the drag force is linear with respect 
to the local differential velocity between the grain and the gas. The 
relation transitions to a power-law regime (Fj^g oc Ai'" "^Av) at inter- 
mediate Reynolds numbers (l < R^ < 800) and becomes quadratic 
at large Reynolds numbers (Ri > 800). When the local concentra- 
tion of dust grains becomes very large (i.e., average distance be- 
tween the particles comparable to the grain size), the coefficient Cd 
should also depend on the local concentration of particles. How- 
ever, this extreme situation is not encountered in astrophysical sit- 
uations. 

Assuming gas molecules interact as hard spheres, the dynamic 
viscosity of the gas can be computed according to ( [Chapman &| 
[Cowling|1970t : 
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where m = Im^i and a-, is the geometric cross section of the 
molecule {as - 2.367 x 10"'^ cm- for H2). The gas mean free 
path /Ig and the kinematic viscosity v of the gas are deduced from 
using 
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3 ASYTROPHYSICAL DUST AND GAS MIXTURES IN 
SPH 

3.1 SPH evolution equations 

The SPH version of the continuity equations Eqs.[T]-[2]are given 
by the density summations for both the gas and the dust phase, 
computed according to: 
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where as in Paper I, the indices a, h, c refer to quantities computed 
on gas particles and i, j, k refer to quantities computed on dust par- 
ticles. The volume filling fraction 9, is defined on a gas particle, a, 
according to 
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where pd „ is the density of dust at the gas particle location, calcu- 
lated using 



Pd.a = 2 injW„j{ha), 



(26) 



where h„ is the smoothing length of the gas particle computed using 
gas neighbours. The local density of dust at the gas location can 
thus be zero (giving 8 = 1) if no dust particles are found within the 
kernel radius computed with the gas smoothing length. Importantly, 
as p and h are mutually dependent, they have to be simultaneously 
calculated for each type of particle, e.g. by the iterative procedure 
described in [Price & Monaghan[ ( [2007] l. 

The SPH equations of motion for the gas and the dust particles, 
corresponding to the SPH translation of Eqs.[3]and|4] are given by 
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for an SPH dust particle. Q. is the usual variable smoothing length 
term 
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neighbours according to: 



£l,ij, - 1 ■ 



6 is defined according to 
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At this stage, no assumptions are made with respect to the func- 
tional form of the drag coefficient K. The evolution of the internal 
energy for an SPH gas particle is given by 



-J- = 7^-^ > , mt (V„ - Vi) • VaWalAK) 



(32) 
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In Paper I, we showed that the total linear and angular momentum 
as well as the total energy are exactly conserved. Thermal coupling 
terms have been neglected in this paper 

3.2 Kernel functions 

Two different kernels are employed to perform the SPH interpola- 
tions. First, a standard bell-shaped kernel W: 

W {r,h) = (q) , (33) 

where h denotes the smoothing lengths of each phases, v the num- 
ber of spatial dimensions and q = \r - r'\/h is the dimensionless 
variable used to calculate the densities and the buoyancy terms. The 
function / is usually the M4 cubic spline kernel ( |Monaghan|2005[ l. 
The drag interpolation is performed using a second kernel D. As 
shown in Paper I, double-hump shaped kernels given by 



D (r, h) ■■ 



(34) 



significantly improve the accuracy of the drag interpolation — for 
the same computational cost — compared to bell-shaped kernels. 
The normalisation constants a for various double hump kernels are 
given in Paper I. We adopt the double hump cubic for the drag terms 
in this paper. 

3.3 Astrophysical drag regimes in SPH 

3.3. 1 Gas viscosity and mean free path 

The drag coefficients Kak involved in Eqs.[27]-|28]and[32]are com- 
puted independently for each pair of any gas particle a and dust 
particle k. We first use the sound speed Cj ii to estimate the viscosity 
Ha on the gas particle a using (see Eq.|20^ 

^ (35) 



640", \ y 

The mean free path is then computed according to Eq.|2T] giving 
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Finally, A„ „ is compared to the quantity 4sk/9 — s^ being the grain 
size of the dust particle — to determine whether the drag coefficient 
of the SPH pair K^k is calculated using the Epstein or the Stokes 
drag regimes. 



3.3.2 Epstein regime 

If '4 Ski 9 < A^a, the drag coefficient K^k is calculated using the Ep- 
stein prescription. Introducing the SPH quantity tf/ak calculated on 
a pair of gas and dust SPH particles and defined by 



(37) 
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Eq.[TT]can be straightforwardly translated to get the drag coefficient 
K„k involved in the SPH drag force 
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where s is the grain radius, is the grain mass and 7 is the adia- 
batic index. Eq.[38]is computationally expensive as it involves ex- 
ponential and error functions. The SPH equivalent of Eq[T7]is given 
by 
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Both Eqs. [38] and [39] reduce to the linear Epstein regime at low 
Mach numbers (equivalent of Eq.|14[( for which the coefficient Kak 
is 
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and to the quadratic drag regime at high Mach numbers (equivalent 
of Eq.[T6l, for which the coefficient Kak is 



Kak = r^P^S- \\ak\- 

mi 



(41) 



3.3.3 Stokes regime 

If 4ij./9 > /ig_o, the drag coefficient Kak is calculated using the 
Stokes prescription (see Eqs. [T8]jl9[ ). The local Reynolds number 
^d.u* is computed for each pair of gas and dust particles using 
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such that the drag coefficient Kak can be computed according to 
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These expressions have been used by [Aylilfe et al.jpOl l\ to com- 
pute the drag on planetesimals in a protoplanetary disc. 



4 TIMESTEPING 

4.1 Explicit timesteping 

The simplest method to evolve the evolution equations for the SPH 
particles is to use an explicit integrator (e.g. the standard Leapfrog). 
The stability of the system is guaranteed provided the timestep re- 
mains smaller than a critical value Ar^. In Paper I, we performed 
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a Von Neumann analysis of the continuous equations, deriving tlie 
explicit time stepping criterion 
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for gas and dust particles, respectively, with the minimum being 
taken over all the particle's neighbours. Although this criterion was 
derived in Paper I for linear drag regimes only, it remains valid even 
for non-linear drag regimes where the drag coefficients depend on 
the differential velocity between the particles, i.e. K^t = Kak (WakD- 



4.2 Implicit timestepping 

When the drag timescale becomes smaller than other time scales in 
the system (e.g. the Courant condition or the orbital timescale), the 
timestep restriction of the explicit methods may become prohibitive 
and implicit methods are required. Monaghan (1997) considered 
the application of two implicit schemes (the first-order Backward- 
Euler and second-order Tischer scheme) to SPH dust-gas mixtures. 
Both schemes are unconditionally stable, but a higher accuracy is 
achieved with second-order schemes. 



4.2.1 Backward-Euler method 

The Backward-Euler scheme applied to the drag interaction be- 
tween SPH dust and gas particles is given by 
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Although the scheme is unconditionally stable, the implicit equa- 
tion with respect to the velocities v"^' must be solved at each time 
step. Direct numerical inversion of this linear system would be pro- 
hibitive given the typical number of neighbour interactions for each 
SPH particle. Thus, approximate or iterative solutions to Eqs.|45|- 
|46|are required. 



4.2.2 Monaghan (1997) scheme 

MonaghMi| jl997| suggested approximating the velocities v"^' of 
Eqs.|45|-|46|using a pairwise treatment in order to preserve the ex- 
act conservation of linear and angular momentum in the SPH for- 
malism. Considering the interaction between the SPH gas particle a 
the dust particle |Monaghan| { 1 997| l introduced pairwise auxiliary 
velocities v defined by: 
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Eqs.|47]and|48]are solved, for a given pair of particles, by taking the 
scalar product by fa, of the difference of the two equations, giving 



l+Af«(m„ + ra,)' 



(49) 



Substituting this expression into Eg. |47| and |48| gives expressions 
for Vo and v,. Iterating this pairwise process by looping over all the 



SPH particles provides an approximate solution for the velocities 
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The main drawback of this method is that the approximation given 
by Eqs.[50]and[5T]is inexact - that is, it provides only an approx- 
imate solution to Eqs.|45]and|46] Furthermore the accuracy of the 
approximation is not known a priori and there is no possibility of 
performing repeated sweeps in order to converge to a more accu- 
rate solution. In practice, we find that the velocities obtained by this 
scheme (for example on the dustybox test) can be significantly in 
error, with no possibility of improving the convergence (for exam- 
ple, by doing several iterations/sweeps). 



4.2.3 Alternative pairwise treatment for linear drag regimes 



We propose a more consistent method for solving Eqs. |45f|46| on a 
given gas or dust particle (a and /, respectively) by sweeping over 
all particle pairs and updating the velocities iteratively according to 
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where v** refers to the improved approximation to v"^' obtained 
after updating each pair and v* to the previous iteration value of v" . 
Eqs.|52]and|53]are solved for each pair of particles by taking the 
dot product of f^, with the difference of the two equations, giving 



(Va/ • f„i)" 



1 + Ar%^ (m„ + m,) 



(54) 



Substituting Eq.|54]in Eqs. [52] and [53] gives the updated velocities 
for the pain Note that during the global sweep over particle pairs 
V* begins as v" at the first iteration but is updated as soon as new 
values become available. 

This pairwise correction ensures that i) both the linear and the 
angular momentum are exactly conserved and ii) the velocities con- 
verge to the correct solution of the implicit scheme given by Eqs.|45| 
and [46] since the last term of Eqs. |52] and |53] tends to zero as the 
number of iterations increases. We thus refine our approximation 
to the solution by performing as many successive iterations as are 
required to reach a suitable convergence criterion. 

4.2.4 Convergence criterion 

We consider that the approximation we obtain from the implicit 
scheme described above is accurate enough when 



(55) 



is satisfied for each particle. Typically, we adopt e = 10""*, which 
ensures that the approximation we make on the time stepping is 
negligible compared to the 0{h^) truncation error of the underlying 
SPH scheme. 
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4.2.5 Implementation into Leapfrog 



4.2.6 Generalisation to non-linear drag regimes 



The Leapfrog scheme is well suited to the evolution of particle 
methods because, for position-dependant forces, it preserves geo- 
metric properties of particle orbits and requires only one evaluation 
per timestep to give second order accuracy. In the standard formu- 
lation, the evolution is computed according to 
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corresponding to Kick, Drift and Kick steps respectively. Adapting 
Leapfrog to deal with velocity dependent forces (e.g. drag) is a pri- 
ori more difficult since for velocity-dependent forces, the last Kick 
is implicit in v' . For our present purposes, this does not present a 
major problem since the drag is already computed implicitly. Split- 
ting the forces into position-dependent (fspn) and drag (fdmg) con- 
tributions, the scheme becomes 
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where the Drag steps represent the implicit updates computed as 
described in Sec. |4.2.3l The disadvantage of Eq.[57]is that two drag 
force evaluations are required, removing one of the advantages of 
the Leapfrog integrator. Inspection of|57|reveals that an alternative 
version that requires only one Drag step can be constructed accord- 
ing to 
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(58) 

where we have combined the drag steps by predicting the velocity 
v^'^. Note that strictly, the Drag step in this method is semi-implicit 
since the force is evaluated using x' rather than x^^'^. However, we 
expect this approximation to be reasonable as at high drag (for 
which the implicit method is designed), the drag mainly changes 
the differential velocity between the fluids and has less of an effect 
on the positions. Care is also required when the timestep changes 
between the steps. We have indicated the correct procedure by spec- 
ifying Afo and Af 1 where Af 1 is the timestep computed based on x' . 

\ requires that f is known at the beginning of the in- 
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Finally, Eq. 

tegration. This can be easily achieved by performing the Drag step 
in Eq. 58 with v''^ = v", A/q = and Afj equal to the timestep 



To extend this alternative pairwise treatment to any non-linear drag 
regime, two additional points have to be considered. Firstly, al- 
though in principle six quantities (v,. for each particle) have to be 
determined for each pair, this can be reduced to a single unknown 
quantity since the drag coefficient depends only on the modulus of 
the differential velocity and the exchange of momentum is directed 
along the line of sight joining the particles. The system of equations 
for a single pair thus reduces to 



where 
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PaPi 



■ + V2. 



(59) 



(60) 



(61) 



calculated using the initial particle positions. 



Secondly, taking the dot product of f„, with the difference of the 
two equations Eqns.|59]and|60]does not lead in general to an equa- 
tion which can be solved analytically. The values of Va, • must 
therefore be determined using a numerical rootfinding procedure 
(we use a Newton-Raphson scheme) before being substituted in 
Eqns.[59]and[60]to determine the velocities for both the gas and the 
dust particles. 

4.2.7 Performance of the implicit scheme 

The computational cost of a timestep with the implicit pairwise 
treatment is more expensive than an explicit timestep since at least 
two iterations have to be performed to ensure that the scheme is 
converged. However, the implicit pairwise treatment will be more 
efficient provided that the number of iterations is much smaller than 
the number of explicit timesteps that would otherwise be required. 

We find in practice that the efficiency of the pairwise treatment 
is mainly determined by the number of iterations required to sat- 
isfy Eq.[55](this aspect was not addressed in the Monaghan](T997] l 
scheme where only one iteration is ever taken in the hope that the 
approximation is sufficiently accurate). The rapidity of the conver- 
gence depends primarily on the ratio r = At/t^ of the timestep 
over the drag stopping time (defined in Eq. (96) of Paper I) and 
on the value of e. Empirically, we have found that, for e = 10"* 
and 1 < r < 10, the implicit pairwise treatment converges effi- 
ciently, the ratio |v''^' - v''|/(minCs) decreasing by ~ two orders 
of magnitude at each iterations. Thus, the implicit pairwise treat- 
ment improves the computational time by a factor of ~ 1-10. How- 
ever, this rapidity of convergence decreases as r increases. At very 
high drag (r > 1000), we find that the implicit scheme becomes 
less efficient than explicit timestepping due to the large number 
of iterations required. A similar behaviour has been found using 
the Gauss-Seidel iterative scheme developed by |Whitehouse et aT] 
l |2005| l to treat SPH radiative transfer in the flux-limited diffusion 
approximation (Bate 2011, private communication), so this issue is 
not specific to the pairwise treatment. 

It is important to note that the computational gain obtained 
with the pairwise scheme does not solve the resolution issue at high 
drag extensively discussed in Paper I. Both of these problems sug- 
gest that a more efficient method for handling high drag regimes is 
required. Such a method is beyond the scope of the present paper. 
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4.2.8 Higher order implicit schemes 

Higher temporal accuracy may be achieved by using second in- 
stead of first order imphcit schemes. The gain in accuracy is ob- 
tained by dividing the drag timestep A? into two half timesteps. 
IMonaghan] jl997[ > suggested the 'Tischer' scheme, where the two 
half timesteps are given by 



Ar/2 



-0.6v Vrat^^ (v"!^- •f„j.)ffltDfl, 

4-^ PaPk \ / 



-0.4v y {y'li^ ■ v„t) r„kD„k, 



(62) 



Ar/2 



and then 



+0 



+0.4v y nil, ^ (v;;,. ■ hi) hiDbi, 



(63) 



1.4v„ - - 0.4v" 



, PaPk 



At/2 



1.4v. - - 0.4v" 



(64) 



= +0.6v y mt (vj+' • ft,) ffaOi, 

V PbPi 



At/2 



(65) 



The last terms of Eqs. |62}{63] correspond to the explicit drag force 
involved in the Forward-Euler scheme. These quantities are com- 
puted form the velocities at the timestep n at the beginning of 
the scheme, as described in Sec. |4.1| Then, the successive deter- 
mination of v"^2 and v" as given by Eqs. |64[ 65 consists of two 
Backward-Euler steps with step size y. They are therefore com- 
puted using our alternative pairwise scheme, until the iterations for 
each half time step have converged. Eqs. |62] and |63] concerns the 
specific case of a linear drag regime, but this scheme can easily be 
extended to non-linear drag regimes as in Sec. |4.2.6] 



5 NUMERICAL TESTS 

5.1 DusTYBOx: Two fluid drag in a periodic box 



The DUSTYBOX problem presented by jLaibe & Price| ( |201 la^ and de- 
scribed in detail in Paper I consists of two fluids in a periodic box 
moving with a differential velocity (Avo = Vrf,o - ^g.o)- This is the 
only test where analytic solutions are known for several functional 
forms corresponding to non-linear drag regimes (see jLaibe & Price] 
|201 la) . These represent the functional forms of the Epstein and 
Stokes prescription. We thus use the dustybox problem to bench- 
mark the accuracy of our algorithm for non-linear drag regimes 
using both explicit and implicit timestepping. 



5.7.7 DUSTYBOX.' jetop 

We set up constant densities pg and and gas pressure in a 3D 
periodic domain x,y,z 6 [0, 1] filled by 20' gas particles set up on 
a regular cubic lattice and 20'' dust particles shifted by half of the 
lattice spacing in each direction (as in Paper I, we verified that the 
results are independent of the offset of the dust lattice). The gas 
sound speed, the gas and the dust densities are set to unity in code 




Figure 1. Dust velocity (solid lines) as a function of time in the dustybox 
test, using 2 X 20' particles, a dust-to-gas ratio of unity and five different 
linear and non-linear drag regimes — quadratic, power-law, linear, third 
order expansion and mixed, from top to bottom — compared to the exact 
solution for each case (long dashed/red lines). The initial velocities are set 
to I'd, I = 1, Vg , = and the time integration is performed using using the 
pairwise implicit treatment described in Sec. [4] The accuracy (< 0.1%) of 
the SPH treatment for dust-gas mixtures is obtained by using the double- 
hump cubic kernel. 



units and no artificial viscosity is applied. The intrinsic dust volume 
is neglected by assuming 9=1. 

Simulations have been performed using both the explicit 
timestepping presented in Paper I and the implicit pairwise 
timestepping described in Sec. |4] For the latter, we verified that 
both the total linear and angular momentum are exactly conserved 
as expected. 



5.7.2 DUSTYBOX.' different drag regimes 

Fig. [T] shows the results of the dustybox test for the five different 
regimes given in Table 1 of |Laibe & Pricelj201 la] l: linear (K = Kg), 
quadratic (A" = TfolAvl), power-law (K = Ko\Av\", with a = 0.4), 
third order expansion (K = Ko[l + ajjAvp], with 03 = 0.5) and 
mixed (K = Kq -y'l 4- aalAvp, with 02 = 5) where we have used 
Kg = I in each case. The analytic solutions are reproduced within 
an accuracy comprised between 0. 1 % and 1 % in every case — both 
linear and non-linear. The implicit scheme was find to converge 
quickly for this problem, requiring no more than two iterations at 
every stage of the evolution in each case. 

The efficiency of the damping for the dustybox problem de- 
creases when the exponent of the drag regime increases, since 
||Av|| < 1. On the contrary, additional non-linear terms give an ad- 
ditional contribution to the drag for the mixed and third order drag 
regimes, leading to a differential velocity that is more efficiently 
damped compared to the linear case. 
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5.2 dustywave: Sound waves in a dust-gas mixture 

The exact solution for linear waves propagating in a dust-gas mix- 
ture (dustywave) was derived in [Laibe & PriceH2011a[ l assuming 
a linear drag regime. Unfortunately, exact solutions prove difficult 
to obtain for the case of non-linear drag. Instead, we have veri- 
fied simply that our simulations of the dustywave problem for non- 
linear drag regimes are converged in both space and time with an 
explicit timestepping scheme. We have then used these results to 
benchmark our simulations using implicit timestepping. 

We run the dustyywave problem for the same five non-linear 
drag regimes used above. Strictly speaking. It should be noted that 
assuming Ko is constant (which we assume for this test problem) 
corresponds to Epstein and Stokes drag only to first order for the 
dustywave problem. 



5.2.1 dustywave.- Setup 

The DUSTYWAVE test is performed in a ID periodic box, placing 
equally spaced particles in the periodic domain x e [0, 1] such that 
the gas and dust densities are unity in code units. We do not apply 
any form of viscosity and the gas sound speed is set to unity. To 
remain in the linear acoustic regime, the relative amplitude of the 
perturbation of both velocity and density are set to lO""*. 



5.2.2 DUSTYWAVE.- Different drag regimes 

Fig. |2] shows the velocity profiles after 5 periods for three drag 
regimes — linear, power law and quadratic — in the ID dusty- 
wave problem using A'o = 1 and a dust-to-gas ratio of unity. The 
solution obtained for the linear drag regime shows an efficiently 
damped perturbation at f = 5, consistent with a stopping time of or- 
der unity. By comparison, the perturbation is only weakly damped 
for the case of the power-law drag regime at the same time. The 
drag is weaker still in the quadratic regime, for which both the dust 
and the gas are mostly decoupled. Indeed, since the drag stopping 
time is a decreasing function of the differential velocity for non- 
linear drag, and the differential velocity is small with respect to the 
sound speed, the damping is inefficient. 

We have also performed dustywave simulations using the third 
order expansion and mixed drag regimes. However, for these cases, 
non-linear terms represent only negligible corrections compared to 
the linear term, thus giving the same results as for the linear case. 



5.3 dustyshock: Two fluid dust-gas sliocks 

The dustyshock problem (Paper I) is a two fluid version of the 
standard |Sod| ([1978} shock tube problem. Results were presented 
in Paper I using a constant drag coefficient K and no heat trans- 
fer between the gas and the dust phase. Here we extend the test to 
non-linear drag regimes. While the evolution during the transient 
stage is dependent on the drag regime, the solution during the sta- 
tionary stage remains unchanged, being a fixed function of the gas 
sound speed and the dust-to-gas ratio. To facilitate the comparison 
between a physical Epstein drag and the case of a constant drag 
coefficient (Paper I), we fix the ratio s^/m^ — s being the grain 
size and the grain mass — to unity in code units, such that the 
Epstein drag coefficient is unity in regions where p„ = I. 



1x10-4 - - 
: t=5 : 

5x10-5 -- - 



-5x10-5 - 
-1x10-4 - 




I \ I I I I I I I \ I 

0.5 1 
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Figure 2. Solution of the ID dustywave problem showing the SPH gas 
(solid black lines) and dust (dashed black) velocities using three different 
drag regimes: linear (top panel), power-law with an exponent of 0.4 (center 
panel) and quadratic (bottom panel). The damping is strongly reduced for 
non-linear drag regimes compared to the linear case. Results are shown af- 
ter 5 wave periods and the linear case (top panel) may be compared to the 
exact analytic solution (red lines). The dust to gas ratio and Kq are set to 1 
in code units. 



5.3.1 DUSTYSHOCK.- Setup 

Equal mass particles are placed in the ID domain x e [-0.5,0.5], 
where for x < we use = = 1, Vg = Vd = and = 1, 
while for X > Opg = = 0.125, Vg = = and P„ = 0.1. We use 
an ideal gas equation of state P = (y - l)pu with y = 5/3. Initial 
particle spacing to the left of the shock in both fluids is Ax = 0.001 
while to the right it is Ax = 0.008, giving 569 equal mass particles 
in each phase. Standard SPH artificial viscosity and conductivity 
terms are applied as in Paper I. 



5.3.2 DUSTYSHOCK.- Different drag regimes 

Fig.|3]illustrates how the transient regime of the dustyshock is af- 
fected when treating the drag with an astrophysical prescription 
where ffie drag coefficient depends on the local density and the 
gas sound speed (right panel) rather than a constant coefficient 
(left panel). Specifically, we use the non-linear Epstein drag regime 
given by Eq. |39] In this dustyshock test, the non-linear terms con- 
stitute a small correction to the linear Epstein drag regime given 
by Eq. [40] Fig. [5] shows that in the case of the Epstein regime, the 
drag is less efficient than for the constant coefficient case, leading 
to a larger (~ by a factor 7) differential velocity between the gas 
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Figure 3. Results of the dustyshock problem with a linear constant drag coefficient (K = I) (left panel) and a non-linear Epstein drag (right panel) where the 
drag coefficients are initially the saine ahead of the shock. The dust-to-gas ratio is set to unity. At f = 0.2, the solutions are in the transient stage where the 
analytic solution is not known. As an indication, the solution for the later stationary stage is shown by the dotted/red lines. The profiles dift'er as the damping 
is less efficient using the non-linear Epstein regime. 



and the dust after t = 0.2 in code units. The dust velocity profile is 
also smoother than in the constant coefficient case. As a result, the 
kinetic energy is less efficiently dissipated by the drag, leading to a 
less sharp peak in the internal energy of the gas. The density profile 
of the dust is also closer to its initial profile behind and ahead of the 
shock. 



5.4 DusTYSEDOv: Two fluid dust-gas blast wave 

The DUSTYSEDOV problem (Paper I) involves the propagation of a 
blast wave in an astrophysical mixture of dust and gas. We adopt 
physical units for this problem, assuming a box size of 1 pc, an am- 
bient sound speed of 2x 10* cm/s and a gas density of po = 6x 10""^^ 
g/cm^ the energy of the blast is 2 x 10" erg and time is measured 
in units of 100 years, roughly corresponding to a supernova blast 
wave propagating into the interstellar medium. We these units, we 
choose the grain size, 0.1/jm and the dust-to-gas ratio, 0.01, to be 
typical of the interstellar medium. In code units, this corresponds to 
an initial drag coefficient of A" = 1 outside the blast radius. As for 
the DUSTYSHOCK, we compare the results using a non-linear Epstein 
drag prescription with with the constant coefficient case described 
in Paper 1. 



5.4.1 DUSTYSEDOV.- Setup 

The problem is set up in a 3D periodic box (x,y,z e [-0.5,0.5]), 
filled by 50'' particles for both the gas and the dust. Gas particles 
are set up on a regular cubic lattice, with the dust particles also on a 
cubic lattice but shifted by half of the lattice step in each direction. 
For shock-capturing, we set q-sph = 1 and jSsph = 2 for the artificial 
viscosity terms and a,, = 1 for the artificial conductivity term. An 
ideal gas equation of state P = iy - i)pu is adopted with y = 5/3. 

The internal energy is distributed of the gas over the particles 
located inside a radius r < rj, where r^, is set to 2h (i.e., the radius of 
the smoothing kernel which for 50^ particles and rj = 1.2 is 0.048). 
In code units the total blast energy is £ = 1, with Pg = 1 and 
Pa = 0.01. For r > rj,, the gas sound speed is set to be 2 x 10"'' in 
code units. 

5.4.2 DUSTYSEDOV.- Different drag regimes 

Figs.[4]and[5]show the evolution of the gas and dust mixture where 
a constant drag coefficient is used (top panels of Fig.|4j compared 
to a drag prescribed by the Epstein regime (bottom panels of Fig.|4] 
Fig. [5|. The gas profiles are similar in both cases since the gas is 
poorly affected by the dust given the low dust-to-gas ratio. How- 
ever, the dust density profiles differ, essentially due to the fact that 
the drag coefficient scales with the sound speed and is thus higher 
in the inner blast region for the Epstein case. Thus, the dust is ef- 
ficiently piled up and accumulates in the gas over-density. As a 
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Figure 5. Cross-section slice showing density in tiie midplane in the 3D dustysedov problem, for both the gas (left panel) and the dust (right panel) at ? = 0. 1 . 
Initially, the dust-to-gas ratio is 0.01 and the drag coefficient is given by the Epstein regime for grains of 0.1/im in size. 50^ SPH particles have been used in 
each phase. 
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Figure 4. Results of the 3D dustysedov test, showing the density in the gas 
(left figure) and dust (right figure) from a Sedov blast wave propagating 
in an astrophysical (1% dust-to-gas ratio) mixture of gas and 0.1/im dust 
grains in a 1 pc box. The drag coefficient is constant (K = 1, top panels) 
or given by the Epstein regime (bottom panels). The low dust-to-gas ratio 
means that the gas is only weakly affected by the drag from the dust, and is 
thus close to the self-similar Sedov solution (dotted/red line). In the Epstein 
case, the drag is much higher inside the blast radius and the dust particles 
are efficiently piled up by ffie passage of the gas over-density. 



result, the dust is cleaned up by the gas in the inner regions of the 
blast, but is more concentrated (by ~ 10%) close to the gas over- 
density than for the constant drag coefficient case. 

The results using either explicit or implicit timestepping were 



found to be indistinguishable. For the Epstein case, we found that 
roughly ten iterations were required for the implicit scheme to con- 
verge on this problem. 

5.5 DUSTYDISC 

The DUSTYDISC problem concerns the evolution of a dusty gas mix- 
ture in a protoplanetary disc (see Paper I for details). For our test 
case, we study how the dust distribution is affected when consider- 
ing a general non-linear Epstein drag instead of the standard linear 
regime. The results obtained when implicitly integrating the non- 
linear drag regime have been found to be similar to benchmark tests 
performed with explicit integration. 



5.5.1 DUSTYDISC- Setup 

We setup 10' gas particles and 10' dust particles in a O.OIMq gas 
disc (with O.OOOIMq of dust) surrounding a IMq star The disc ex- 
tends from 10 to 400 AU. Both gas and dust particles are placed 
using a Monte-Carlo setup such that the surface density profiles of 
both phases are X (r) oc r"' . The radial profile of the gas tempera- 
ture is taken to be T {r) oc r""-^ with a flaring H/r = 0.05 at 100 
AU. One code unit of time corresponds to 10^ yrs. 



5.5.2 DUSTYDISC Evolution of the particles 

Fig.|6]shows a face-on view of a protoplanetary disc, integrating the 
linear Epstein regime (left panel) and the full non-linear Epstein 
drag (right panel). The dust distributions are not found to exhibit 
significant discrepancies. In the non-linear drag regime case how- 
ever, the dust distribution is slightly smoother since the drag (and 
thus, the coupling with the gas phase) is more efficient. 

Fig.|7]compares the vertical motion of a dust grain initially lo- 
cated at z = zo using the linear (explicit integration) and the full non 
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Figure 6. Rendering of the density for tlie dust of a typical T-Tauri Star protoplanerary disc using 2 X 10^ SPH particles, using an explicit time integration in 
the linear Epstein regime (left panel) and an implicit integrator in the full non-linear Epstein drag regimes (right panel). 

with respect to the differential velocity between the gas and the 
dust. 

Particular attention has been paid to developing an implicit 
timestepping scheme to efficiently simulate the case of high drag, 
extending the scheme proposed by |Monaghan| (|T997j which we 
found to be unsatisfactory. We have presented a new pairwise im- 
plicit scheme that, like the |Monaghan| ( (T997] l scheme, preserves the 
exact conservation of linear and angular momentum but, unlike the 
[Monaghan^ ( , 1 997 1 scheme, i) provides control over the accuracy of 
the iterative procedure and ii) can incorporate non-linear terms for 
both Epstein and Stokes drag. We found that when the ratio r be- 
tween the the timestep and the drag stopping time is 1 < r < 1000, 
the implicit timestepping is faster than a standard explicit integra- 
tion. However, at higher values of r, the algorithm is less efficient. 

The accuracy of the generalised algorithm is benchmarked 
against the suite of test problems presented in Paper I. In particular, 
the solutions obtained for the dustybox problem are compared to 
their known analytic solutions for a large range of non-linear drag 
regimes and the solutions of the dustywave, dustyshock, dustyse- 
Dov and dustydisc problems are benchmarked against converged 
results obtained with explicit timestepping. 

The two key issues addressed in this paper complete the study 
of our algorithm developed in Paper I. Our intention is to apply it 
to various astrophysical problems involving gas and dust mixtures 
in star and planet formation. A first application is given in .Ayliffe] 
|etaL] l |20n) . 
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Figure 7. Vertical settling of a dust grain (1cm in size) initially located at 
ro = 100 AU and zo = 2 AU (solid/black), integrating implicitly the non- 
linear Epstein regime. SPH results are compared to the explicit integration 
of the linear Epstein regime (dashed/red) and the estimation given by the 
damped harmonic oscillator approximation (pointed/red). In the full non- 
linear drag regime, the settling is more efBcient than for the linear case 
since the vertical oscillations in the dust motion reaches a fraction zo/H of 
the sound speed. 

linear (implicit integration) Epstein regimes. In the full non-linear 
case, the settling is more efficient since the vertical differential ve- 
locity between the dust grains and the gas in the mid plane of the 
disc reaches a fraction zo/H of the sound speed, meaning that the 
non-linear terms are no longer negligible. 



6 CONCLUSIONS 

We have extended the SPH formalism for two-ffuid dust and gas 
mixtures developed in Paper I to handle the drag regimes usually 
encountered in a large range of astrophysical contexts. Specifically, 
our algorithm is now designed to treat the dynamics of grains sur- 
rounded by a dilute medium (Epstein regime) or dense fluid (Stokes 
regime), for which the drag force can be either linear or non-linear 
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